R für die Sozialwissenschaften - Teil 2

Jan-Philipp Kolb

23 Juni 2017

ggplot und ggmap

Das Paket ggplot2

Basiseinführung ggplot2

<www.r-bloggers.com/basic-introduction-to-ggplot2/>

install.packages("ggplot2")
library(ggplot2)

Der diamonds Datensatz

head(diamonds)
carat cut color clarity depth table price x y z
0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43
0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31
0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31
0.29 Premium I VS2 62.4 58 334 4.20 4.23 2.63
0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75
0.24 Very Good J VVS2 62.8 57 336 3.94 3.96 2.48

Wie nutzt man qplot

# histogram
qplot(depth, data=diamonds)

Ein Balkendiagramm

qplot(cut, depth, data=diamonds)

Ein weiteres Balkendiagramm

qplot(factor(cyl), data=mtcars,geom="bar")

Boxplot

qplot(data=diamonds,x=cut,y=depth,geom="boxplot")

Scatterplot

# scatterplot
qplot(carat, depth, data=diamonds)

Farbe hinzu:

qplot(carat, depth, data=diamonds,color=cut)

Trendlinie hinzufügen

myGG<-qplot(data=diamonds,x=carat,y=depth,color=carat) 
myGG + stat_smooth(method="lm")

Graphik drehen

qplot(factor(cyl), data=mtcars, geom="bar") + 
coord_flip()

Wie nutzt man ggplot

ggplot(diamonds, aes(clarity, fill=cut)) + geom_bar()

Farben selber wählen

Es wird das Paket RColorBrewer verwendet um die Farbpalette zu ändern

install.packages("RColorBrewer")
library(RColorBrewer)
myColors <- brewer.pal(5,"Accent")
names(myColors) <- levels(diamonds$cut)
colScale <- scale_colour_manual(name = "cut",
                                values = myColors)

http://stackoverflow.com/questions/6919025/

Eine Graphik mit den gewählten Farben

p <- ggplot(diamonds,aes(carat, depth,colour = cut)) + 
  geom_point()
p + colScale

Speichern mit ggsave

ggsave("Graphik.jpg")

Gliederung

Arten von räumlichen Daten:

Das R-paket ggmap wird im folgenden genutzt um verschiedene Kartentypen darzustellen.

Mit qmap kann man eine schnelle Karte erzeugen.

Straßenkarten

Installieren des Paketes

devtools::install_github("dkahle/ggmap")
devtools::install_github("hadley/ggplot2")
install.packages("ggmap")

Paket ggmap - Hallo Welt

library(ggmap)

Und schon kann die erste Karte erstellt werden:

qmap("Mannheim")

Karte für eine Sehenswürdigkeit

BBT <- qmap("Berlin Brandenburger Tor")
BBT

Karte für einen ganzen Staat

qmap("Germany")

Ein anderes zoom level

qmap("Germany", zoom = 6)

Hilfe bekommen wir mit dem Fragezeichen

?qmap

Verschiedene Abschnitte in der Hilfe:

Die Beispiele in der Hilfe

Ausschnitt aus der Hilfe Seite zum Befehl qmap:

qmap Example

qmap Example

Das Beispiel kann man direkt in die Konsole kopieren:

# qmap("baylor university")
qmap("baylor university", zoom = 14)
# und so weiter

Ein anderes zoom level

qmap("Mannheim", zoom = 12)

Näher rankommen

qmap('Mannheim', zoom = 13)

Ganz nah dran

qmap('Mannheim', zoom = 20)

ggmap - maptype satellite

qmap('Mannheim', zoom = 14, maptype="satellite")

ggmap - maptype satellite zoom 20

qmap('Mannheim', zoom = 20, maptype="hybrid")

ggmap - maptype hybrid

qmap("Mannheim", zoom = 14, maptype="hybrid")

Terrain/physical maps

ggmap - terrain map

qmap('Schriesheim', zoom = 14,maptype="terrain")

Abstrahierte Karten (http://www.designfaves.com)

New York

New York

ggmap - maptype watercolor

qmap('Mannheim', zoom = 14,maptype="watercolor",source="stamen")

ggmap - source stamen

qmap('Mannheim', zoom = 14,
 maptype="toner",source="stamen")

ggmap - maptype toner-lite

qmap('Mannheim', zoom = 14,
 maptype="toner-lite",source="stamen")

ggmap - maptype toner-hybrid

qmap('Mannheim', zoom = 14,
 maptype="toner-hybrid",source="stamen")

ggmap - maptype terrain-lines

qmap('Mannheim', zoom = 14,
 maptype="terrain-lines",source="stamen")

Graphiken speichern

RstudioExport

RstudioExport

ggmap - ein Objekt erzeugen

MA_map <- qmap('Mannheim', 
               zoom = 14,
               maptype="toner",
               source="stamen")

Geokodierung

Geocoding (…) uses a description of a location, most typically a postal address or place name, to find geographic coordinates from spatial reference data …

Wikipedia - Geocoding

library(ggmap)
geocode("Mannheim",source="google")
lon lat
8.463182 49.48608

Latitude und Longitude

LatLon

LatLon

http://modernsurvivalblog.com

Koordinaten verschiedener Orte in Deutschland

cities lon lat
Hamburg 9.993682 53.55108
Koeln 6.960279 50.93753
Dresden 13.737262 51.05041
Muenchen 11.581981 48.13513

Reverse Geokodierung

Reverse geocoding is the process of back (reverse) coding of a point location (latitude, longitude) to a readable address or place name. This permits the identification of nearby street addresses, places, and/or areal subdivisions such as neighbourhoods, county, state, or country.

Quelle: Wikipedia

revgeocode(c(48,8))
## [1] "Unnamed Road, Somalia"

Die Distanz zwischen zwei Punkten

mapdist("Q1, 4 Mannheim","B2, 1 Mannheim")
##             from             to   m    km     miles seconds  minutes
## 1 Q1, 4 Mannheim B2, 1 Mannheim 749 0.749 0.4654286     215 3.583333
##        hours
## 1 0.05972222
mapdist("Q1, 4 Mannheim","B2, 1 Mannheim",mode="walking")
##             from             to   m    km     miles seconds minutes  hours
## 1 Q1, 4 Mannheim B2, 1 Mannheim 546 0.546 0.3392844     423    7.05 0.1175

Eine andere Distanz bekommen

mapdist("Q1, 4 Mannheim","B2, 1 Mannheim",mode="bicycling")
##             from             to   m    km    miles seconds  minutes
## 1 Q1, 4 Mannheim B2, 1 Mannheim 555 0.555 0.344877     215 3.583333
##        hours
## 1 0.05972222

Geokodierung - verschiedene Punkte von Interesse

POI1 <- geocode("B2, 1 Mannheim",source="google")
POI2 <- geocode("Hbf Mannheim",source="google")
POI3 <- geocode("Mannheim, Friedrichsplatz",source="google")
ListPOI <-rbind(POI1,POI2,POI3)
POI1;POI2;POI3
##        lon      lat
## 1 8.462844 49.48569
##        lon      lat
## 1 8.469879 49.47972
##        lon      lat
## 1 8.475208 49.48326

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),
data = ListPOI)

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),col="red",
data = ListPOI)

ggmap - verschiedene Farben

ListPOI$color <- c("A","B","C")
MA_map +
geom_point(aes(x = lon, y = lat,col=color),
data = ListPOI)

ggmap - größere Punkte

ListPOI$size <- c(10,20,30)
MA_map +
geom_point(aes(x = lon, y = lat,col=color,size=size),
data = ListPOI)

Eine Route von Google maps bekommen

from <- "Mannheim Hbf"
to <- "Mannheim B2 , 1"
route_df <- route(from, to, structure = "route")

Mehr Information

http://rpackages.ianhowson.com/cran/ggmap/man/route.html

Eine Karte mit dieser Information zeichnen

qmap("Mannheim Hbf", zoom = 14) +
  geom_path(
    aes(x = lon, y = lat),  colour = "red", size = 1.5,
    data = route_df, lineend = "round"
  )

Wie fügt man Punkte hinzu

http://i.stack.imgur.com

Cheatsheet

https://www.rstudio.com/

Cheatsheet

Cheatsheet

Resourcen und Literatur

Take Home Message

Was klar sein sollte:

Die lineare Regression

Die lineare Regression

Maindonald - Data Analysis

Lineare Regression in R - Beispieldatensatz

John H. Maindonald and W. John Braun

DAAG - Data Analysis and Graphics Data and Functions

install.packages("DAAG")
library("DAAG")
data(roller)

help on roller data:

?roller

Das lineare Regressionsmodell in R

Schätzen eines Regressionsmodells:

roller.lm <- lm(depression ~ weight, data = roller)

So bekommt man die Schätzwerte:

summary(roller.lm)
## 
## Call:
## lm(formula = depression ~ weight, data = roller)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.180 -5.580 -1.346  5.920  8.020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.0871     4.7543  -0.439  0.67227   
## weight        2.6667     0.7002   3.808  0.00518 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.735 on 8 degrees of freedom
## Multiple R-squared:  0.6445, Adjusted R-squared:  0.6001 
## F-statistic:  14.5 on 1 and 8 DF,  p-value: 0.005175

Falls das Modell ohne Intercept gesch?tzt werden soll:

lm(depression ~ -1 + weight, data = roller)
## 
## Call:
## lm(formula = depression ~ -1 + weight, data = roller)
## 
## Coefficients:
## weight  
##  2.392

Summary des Modells

summary(roller.lm)
## 
## Call:
## lm(formula = depression ~ weight, data = roller)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.180 -5.580 -1.346  5.920  8.020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.0871     4.7543  -0.439  0.67227   
## weight        2.6667     0.7002   3.808  0.00518 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.735 on 8 degrees of freedom
## Multiple R-squared:  0.6445, Adjusted R-squared:  0.6001 
## F-statistic:  14.5 on 1 and 8 DF,  p-value: 0.005175

R arbeitet mit Objekten

predict(roller.lm) # Vorhersage
##         1         2         3         4         5         6         7 
##  2.979669  6.179765  6.713114 10.713233 12.046606 14.180002 14.980026 
##         8         9        10 
## 18.180121 24.046962 30.980502
resid(roller.lm) # Residuen
##          1          2          3          4          5          6 
## -0.9796695 -5.1797646 -1.7131138 -5.7132327  7.9533944  5.8199976 
##          7          8          9         10 
##  8.0199738 -8.1801213  5.9530377 -5.9805017

Residuenplot

plot(roller.lm,1)

Residuenplot

plot(roller.lm,2)

Regressionsdiagnostik mit Basis-R

Ein einfaches Modell

N <- 5
x1 <- rnorm(N)
y <- runif(N)
par(mfrow=c(1,2))
plot(density(x1))
plot(density(y))

Modellvorhersage machen

mod1 <- lm(y~x1)
pre <- predict(mod1)
y
## [1] 0.09257606 0.72996711 0.38776661 0.43490811 0.88667362
pre
##         1         2         3         4         5 
## 0.5034355 0.4927091 0.5027524 0.5193458 0.5136487

Regressionsdiagnostik mit Basis-R

plot(x1,y)
abline(mod1)
segments(x1, y, x1, pre, col="red")

Beispieldaten Luftqualität

library(datasets)
?airquality

Das visreg-Paket

Ein Modell wird auf dem airquality Datensatz geschätzt

install.packages("visreg")
library(visreg)
fit <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality)
summary(fit)
## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind + Temp, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.485 -14.219  -3.551  10.097  95.619 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
## Solar.R       0.05982    0.02319   2.580  0.01124 *  
## Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
## Temp          1.65209    0.25353   6.516 2.42e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.18 on 107 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.6059, Adjusted R-squared:  0.5948 
## F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

Visualisierung

visreg(fit)

Und dann mit visreg visualisiert.

visreg(fit, "Wind", type = "contrast")

Visualisierung mit dem Paket visreg

visreg(fit, "Wind", type = "contrast")

Das visreg-Paket

visreg(fit, "Wind", type = "conditional")

Regression mit Faktoren

Mit visreg können die Effekte bei Faktoren visualisiert werden.

airquality$Heat <- cut(airquality$Temp, 3, 
    labels=c("Cool", "Mild", "Hot"))
fit.heat <- lm(Ozone ~ Solar.R + Wind + Heat, 
    data = airquality)
summary(fit.heat)
## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind + Heat, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.473 -12.794  -2.686   8.461 107.035 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 48.27682    8.83072   5.467 3.07e-07 ***
## Solar.R      0.06524    0.02145   3.042  0.00296 ** 
## Wind        -3.49803    0.59042  -5.925 3.94e-08 ***
## HeatMild     9.05367    4.88257   1.854  0.06648 .  
## HeatHot     43.13970    5.98618   7.207 8.79e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.9 on 106 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.6554, Adjusted R-squared:  0.6424 
## F-statistic:  50.4 on 4 and 106 DF,  p-value: < 2.2e-16

Effekte von Faktoren

par(mfrow=c(1,2))
visreg(fit.heat, "Heat", type = "contrast")
visreg(fit.heat, "Heat", type = "conditional")

Das Paket visreg - Interaktionen

airquality$Heat <- cut(airquality$Temp, 3, 
labels=c("Cool", "Mild", "Hot"))
fit <- lm(Ozone ~ Solar.R + Wind * Heat, data = airquality)
summary(fit)
## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind * Heat, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.472 -11.640  -1.919   7.403 102.428 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.48042   17.38219   0.258 0.797102    
## Solar.R        0.07634    0.02137   3.572 0.000538 ***
## Wind           0.05854    1.34860   0.043 0.965458    
## HeatMild      56.72928   18.53110   3.061 0.002805 ** 
## HeatHot       94.68619   18.61619   5.086 1.63e-06 ***
## Wind:HeatMild -4.11933    1.57108  -2.622 0.010054 *  
## Wind:HeatHot  -4.88125    1.74358  -2.800 0.006101 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.28 on 104 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.6825, Adjusted R-squared:  0.6642 
## F-statistic: 37.26 on 6 and 104 DF,  p-value: < 2.2e-16

Steuern der Graphikausgabe mittels layout

visreg(fit, "Wind", by = "Heat",layout=c(3,1))

Das Paket visreg - Interaktionen overlay

fit <- lm(Ozone ~ Solar.R + Wind * Heat, data = airquality)
visreg(fit, "Wind", by="Heat", overlay=TRUE, partial=FALSE)

Das Paket visreg - visreg2d

fit2 <- lm(Ozone ~ Solar.R + Wind * Temp, data = airquality)
visreg2d(fit2, "Wind", "Temp", plot.type = "image")

Das Paket visreg - surface

visreg2d(fit2, "Wind", "Temp", plot.type = "persp")

Linkliste - lineare Regression

Die logistische Regression

Agresti - Categorical Data Analysis (2002)

Faraway Bücher zu Regression in R

Binäre AVs mit glm

Beispieldaten für die logistische Regression

install.packages("HSAUR")
library("HSAUR")
data("plasma", package = "HSAUR")

Logistische Regression mit R

cdplot(ESR ~ fibrinogen, data = plasma)

Logistische Regression mit R

plasma_glm_1 <- glm(ESR ~ fibrinogen, data = plasma, 
                    family = binomial())

Weitere Beispieldaten

install.packages("faraway")
library("faraway")
data(orings)
temp damage
53 5
57 1
58 1

Generalisierte Regression mit R - weitere Funktionen

probitmod <- glm(cbind(damage,6-damage) ~ temp, 
    family=binomial(link=probit), orings)
modp <- glm(Species ~ .,family=poisson,gala)
library("MASS")
house.plr<-polr(Sat~Infl,weights=Freq,data=housing)

Linkliste - logistische Regression

Mehrebenenmodelle

Wie sehen die Daten aus?

Andres Gutierrez - Multilevel Modeling of Educational Data using R

Beispiel Mehrebenenmodelle

Untersuchungsgegenstand

Fragen hierzu

Beispiel in Goldstein (2010), Kapitel 1.2

Evaluierung der Effektivität von Schulen

Mehrebenen-Modelle:

Unterscheidung

Bibliotheken

install.packages("lme4")
install.packages("sjPlot")
library(ggplot2)
library(gridExtra)
library(lme4)
library(sjPlot)
library(dplyr)

Beispieldaten

mlexdat <- read.csv("https://github.com/Japhilko/RSocialScience/raw/master/data/mlexdat.csv") 
X SES Score ID
1 18.62733 -55.120574 A
2 33.64915 -92.375273 A
3 22.26931 -48.783447 A
4 36.49052 38.099329 A
5 38.21402 339.701754 A
6 11.36669 2.286978 A

Formalistisch

$$ y_{ij} = {j} + {ij}

_{j} = 0 + u{j} $$ Die Gesamtvariation wird in zwei Teile zerlegt:

Der R-code für dieses Nullmodell

HLM0 <- lmer(Score ~ (1 | ID), data = mlexdat)

Nullmodell Ergebnis

coef(HLM0)
## $ID
##   (Intercept)
## A     45.7893
## B    430.7218
## C   1182.1760
## D   2145.2329
## E   3489.1408
## 
## attr(,"class")
## [1] "coef.mer"
summary(HLM0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ (1 | ID)
##    Data: mlexdat
## 
## REML criterion at convergence: 7130.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74559 -0.69317 -0.00757  0.68337  2.96511 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 1931758  1389.9  
##  Residual               87346   295.5  
## Number of obs: 500, groups:  ID, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   1458.6      621.7   2.346

Interpratation des Nullmodells

100 * 87346 / (87346 + 1931757)
## [1] 4.32598

Ein weiteres Modell

\[y_{ij} = \alpha_{j} + \beta_{j} * SES_{ij} + \varepsilon_{ij}\]

\[\alpha_{j} = \alpha_0 + u_{j}\] \[\beta{j} = \beta_0 + v_{j}\]

Rcode für dieses Modell

HLM1 <- lmer(Score ~ SES + (SES | ID), data = mlexdat)
coef(HLM1)
## $ID
##   (Intercept)        SES
## A    36.46401  0.3798185
## B    37.21549  9.7596237
## C    38.10719 20.8897245
## D    38.85566 30.2320132
## E    39.70159 40.7907386
## 
## attr(,"class")
## [1] "coef.mer"
summary(HLM1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ SES + (SES | ID)
##    Data: mlexdat
## 
## REML criterion at convergence: 6742.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.83274 -0.64740  0.02662  0.69063  2.67309 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  ID       (Intercept)     1.65   1.285      
##           SES           257.09  16.034  1.00
##  Residual             40400.24 200.998      
## Number of obs: 500, groups:  ID, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   38.069     45.863   0.830
## SES           20.410      7.236   2.821
## 
## Correlation of Fixed Effects:
##     (Intr)
## SES -0.119

# 1% - BS variance
# 99% - WS variance
100 * 40400.24 / (40400.24 + 257.09 + 1.65)
## [1] 99.36363
# Percentage of variation explained by SES between schools
1 - ((257.09 + 1.65) / 1931757)
## [1] 0.9998661
# Percentage of variation explained by SES within schools
1 - (40400.24 / 87346)
## [1] 0.5374689

Auf der einen Seite stellen wir fest, dass die SES 99 Prozent der Unterschiede zwischen den Schulen erklärt; Auf der anderen Seite erklärt die SES 53 Prozent der Abweichungen innerhalb der Schulen.

Was heißt das? Schulsegregation

Das heißt, wohlhabende Studenten gehören zu reichen Schulen, und arme Studenten gehören zu armen Schulen.

Darüber hinaus übertreffen wohlhabende Studenten weit über die Leistung der armen Studenten.

Multilevel Model Specification

library(lme4)
library(mlmRev)
names(Exam)
##  [1] "school"   "normexam" "schgend"  "schavg"   "vr"       "intake"  
##  [7] "standLRT" "sex"      "type"     "student"

random intercept, fixed predictor in individual level

lmer(normexam ~ standLRT + (1 | school), data=Exam)

random intercept, random slope

Das nächste Modell, das angegeben wird, ist ein Modell mit einem zufälligen Intercept auf individueller Ebene und ein Prädiktor, der zwischen Gruppen variieren darf. Mit anderen Worten, die Wirkung der Hausaufgaben auf die Partitur auf einem Mathe-Test variiert zwischen den Schulen. Um dieses Modell zu schätzen, wird das ‘1’, das den Intercept im zufälligen Teil der Modellspezifikation angibt, durch die Variable ersetzt, von der wir den Effekt zwischen den Gruppen variieren wollen.

Varying intercept model

MLexamp.6 <- lmer(extro ~ open + agree + social + (1 | school), data = lmm.data)

Varying slope model

MLexamp.9 <- lmer(extro ~ open + agree + social + (1 + open | school/class), data = lmm.data)

Cross-Level-Interaktionseffekt

Paket lmer

lmer(y ~ 1 + (1 | subjects), data=data)
# nlme
lme(y ~ 1, random = ~ 1 | subjects, data=data)